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Because of the potentially large number of important applications of nonlinear optics, researchers 
have expended a great deal of effort to optimize the second-order molecular nonlinear-optical re- 
sponse, called the hyperpolarizability. The focus of our present studies is the intrinsic hyperpolariz- 
ability, which is a scale-invariant quantity that removes the effects of simple scaling, thus being the 
relevant quantity for comparing molecules of varying sizes. Past theoretical studies have focused 
on structural properties that optimize the intrinsic hyperpolarizability, which have characterized 
the structure of the quantum system based on the potential energy function, placement of nuclei, 
geometry, and the effects of external electric and magnetic fields. Those previous studies focused 
on single-electron models under the influence of an average potential. In the present studies, we 
generalize our calculations to two-electron systems and include electron interactions. As with the 
single-electron studies, universal properties are found that are common to all systems - be they 
molecules, nanoparticles, or quantum gases - when the hyperpolarizability is near the fundamental 
limit. 



I. INTRODUCTION 

Materials with large nonlinear-optical suscep- 
tibilities are central for optical applications such 
as telecommunications, [l[ three-dimensional nano- 
photolithography, 0, |j| and making new materials Q 
for novel cancer therapies. [|| Semi-empirical theoretical 
modeling of complex molecules, with feedback from 
experiments, has resulted in a better understanding of 
how to make materials with larger nonlinear-optical 
response, [a [Zj 

Accurate theoretical modeling of a variety of spe- 
cific molecules is a critical part in the development of 
structure-property relationships. [§4ll1| In contrast, our 
work seeks to build a broad theoretical understanding of 
the nonlinear-optical response that applies to all quan- 
tum systems. [12j] of which molecules are a small subset. 
The uniqueness of this approach is that it seeks to iden- 
tify universal properties that are common to all systems 
sharing a nonlinear response that is close to the funda- 
mental limit. (l3| In the present work, we show that in- 
teractions between electrons do not have an effect on the 
universal properties, and that scaling due to the addition 
of an electron is as predicted from quantum limit theory. 

For electric fields that are small compared with in- 
ternal molecular fields, the induced dipole moment of a 
molecule, p, can be expressed as a series in the applied 
electric field, E, 

p = aE + (3E 2 +jE 3 + . . . , (1) 

where a is the polarizability, /3, the hyperpolarizability, 
7 the second hyperpolarizability, and so on. In general, 
these coefficients are tensors that depend on the frequen- 
cies of the applied electric field. However, our studies 
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focus on (1) the non-resonant regime, were all frequen- 
cies are small compared with the Bohr frequencies of the 
molecule; and, (2) the xxx tensor component of /3, where 
x is the axis along which the hyperpolarizability is the 
largest. 

We consider only the electronic hyperpolarizability in 
the off-resonant limit - where the nuclei are assumed not 
to contribute. Chernyak and coworkers showed that in 
second and third harmonic generation experiments, when 
the wavelengths are tuned below the excited electronic 
state energy and above the nuclear excitation energies, 
the purely electronic hyperpolarizabilities account for 90- 
95% of the total. [l4| Our theoretical analysis and ex- 
perimental results come to the same conclusion. [IH, [l6[ 
Bishop and coworkers point out that while this may 
be true for harmonic generation hyperpolarizabilties, vi- 
bronic contributions can be dominant when any of the 
frequencies vanish or in self action effects where nega- 
tive frequencies cancel positive ones. [131 The literature 
is extensive on the topic of conditions when the vibronic 
contributions are dominant. [l8l - [2~fl | 

In the theoretical studies presented here, we focus 
specifically on the zero-frequency limit of the electronic 
response. Our results can be compared with second- 
harmonic generation measurements of the off-resonant 
hyperpolarizability, which is dominated by the electronic 
response, upon extrapolation to zero frequency 22] using 
a two- level dispersion model. [23] All future references to 
hyperpolarizabilities refer to the first electronic hyperpo- 
larizability at zero-frequency. We stress that our goal 
is to understand the zero-frequency electronic response, 
independent of the subtleties associated with measuring 
it. 

The hyperpolarizability can be calculated for any 
quantum system with the dipole energy due to the ap- 
plied electric field as a perturbation. This results in a 
sum-over-states (SOS) expansion, [24] which is expressed 
in terms of the position matrix elements and energy 
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eigenvalues of the unperturbed molecule. 

The fundamental limit of the hyperpolarizability is de- 
fined as the upper bound determined from fundamental 
quantum principles, without approximation. Thus, the 
fundamental limit of the hyperpolarizability is an abso- 
lute quantity that applies to any quantum system that 
is described by the Schrodinger Equation. Rather than 
being a single number, the limit is a function of the en- 
ergy difference between ground and first excited state, 
E w = Ei — E , and the number of electrons. 

This upper bound is calculated using the Thomas- 
Reiche-Kuhn sum rules, [25j which are exact and follow 
directly from the Schrodinger equation. They relate the 
energy eigenvalues and position matrix elements to each 
other, and are used to simplify the SOS expansion of the 
hyperpolarizability. [24] This approach was applied to the 
calculation the fundamental limit of the off-resonance hy- 
perpolarizability, which yields[26] 
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where N is the number of electrons, Eiq the energy differ- 
ence between the first excited state and the ground state, 
Eiq = Ei — Eq, and m is the electron mass. The fun- 
damental limit of the second hyperpolarizability has also 
been calculated. [27f These calculations show that there 
is a large gap between the hyperpolarizabilities and the 
fundamental limit. (26l - [3ll | The existence of such a limit 
is both scientifically interesting and of practical utility 
because it provides a target for making optimized mate- 
rials. 

One approach to the development of structure- 
property relationships is to correlate the nonlinear- 
optical response with variations in the structure. How- 
ever, it is difficult to decouple structural properties from 
simple scaling - i.e. the change in the effective size of a 
molecule with a change in its structure. It can be shown 
that the intrinsic hyperpolarizability pint, which is ob- 
tained by dividing the actual hyperpolarizability /3 by 
the fundamental limit PmAX- 



Pint — P/PMAX > 



(3) 



is invariant to simple scaling. [13| A comparison of the 
intrinsic hyperpolarizabilties of molecules can be used to 
determine the properties that are relevant for optimizing 
the nonlinear-optical response. 

Past studies have considered the optimization of the 
intrinsic hyperpolarizability by varying the shape of 
the potential energy function, [32|, |33[ which lead to the 
new paradigm of large-/? molecules using modulation 
of conju gati on for which there is some experimental 
evidence; [34 135| studying the effect of geometry by vary- 
ing the charges and positions of point nuclei; 36] and the 
effect of applying an arbitrary electromagnetic field. [37j 
All of these studies find that when the hyperpolariz- 
ability of a quantum system is near the fundamental 
limit, certain universal properties result. For example, 



there are many one-dimensional potentials - each of dis- 
tinct shape, but different from each other - whose op- 
timized intrinsic hyperpolarizability is 0.709. For these 
systems, the ratio of second to first excited state energy 
is E = E10/E20 ~ 0.48. Furthermore, typically over 
90% of the optimized hyperpolarizability is due to two 
excited states. Thus, all optimized materials can be well 
approximated by a three-level model. 

When maximizing the off-resonance hyperpolarizabil- 
ity of a one-dimensional molecule, planar molecules made 
of point nuclei, or systems with externally applied elec- 
tric and magnetic fields, one-electron models were unable 
to get intrinsic hyperpolarizabilities that exceed about 
0.708. However, Monte Carlo simulations that randomly 
assign energies and position matrix elements under the 
constraint that the sum rules are obeyed result in j3i n t 
that approaches unity. [38| Since such large values have 
never been observed under any Hamiltonian-based simu- 
lation, this suggests that one may need to invoke exotic 
Hamiltonians to reach the fundamental limit. 

The values of /3;„t « 0.708 that we have obtained in our 
computer simulations are much higher than any that have 
been measured for real molecules. For example, it was 
reported that a new molecule with asymmetric modula- 
tion of conjugation has a measured value of (3i nt = 0.048 
[3^ |. but this is still more than a factor of ten lower 
than what our simulations have produced. It is natu- 
ral to ask whether our simulations give unrealistically 
large (3i nt values simply because we are using one-electron 
models, [ti] 

For a system of two non-interacting electrons, the value 
of (3 is twice that of a system with one electron. But 
Pmax scales as 2 3 / 2 due to the factor TV 3 / 2 in Equation 
Thus the highest value of fii n t that we will be able 
to get with the intrinsic hyperpolarizability of the two 
electron system is about .708/ y/2 ss 0.5. As the num- 
ber of non-interacting electrons is increased, the intrinsic 
hyperpolarizability decreases. For a 10-electron system, 
this would reduce the intrinsic hyperpolarizability by a 
factor of 30, the gap observed between the best molecules 
and the fundamental limit. 

In the present work we consider two-electron models 
and ask whether it is possible to produce values of P int 
that are substantially greater than 0.5 when including 
electron-electron interactions. We will show that fii n t 
of the two-electron system is as large as about 0.708. 
Furthermore, we find that the same universal properties 
hold as in the one-electron case. Thus, the high j3i n t 
values that we obtained in previous work are not simply 
an artifact of one-electron models. More importantly, 
this suggests that one can use simple one-electron models 
when painting our understanding of the nonlinear-optical 
response with the broadest of brushes. 
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II. THEORY 

We use a two-particle Hamiltonian that includes a one- 
dimensional potential energy function, V{x), that acts on 
both electrons, and an interaction term between the elec- 
trons. Note that past studies in two dimensions show the 
same universal behavior as in one dimension, so we find 
no particular utility in studying these other classes. If 
electron interactions do not change the universal proper- 
ties in ID, it is unlikely that they will have any effect in 
2D. 

The potential function is modeled by a piecewise linear 
function that has 19 degrees of freedom in the compu- 
tational region. With 39 degrees of freedom, the results 
are nearly the same; but, the computations are more time 
consuming. We thus use 19 degrees of freedom because it 
provides a good comprmise between accuracy and speed. 

Starting from a given potential function, we use the 
Nelder-Mead simplex algorithm [40] to vary the shape 
of the potential to maximize Pint- Since there are 
19 degrees of freedom, we are maximizing the hyper- 
polarizability over a 19-dimensional space. We have 
three ways of computing (SOS. pi} dipole free.[4ll| and 
non-perturbative[42j), all of which require solving the 
Schrodinger eigenvalue problem for the given potential 
(and in some cases also for neighboring potentials). We 
solve the eigenvalue problem numerically using the finite 
element method. [43| 

Once we have solved the eigenvalue problem, we can 
compute transition moments and then obtain by the 
standard Orr and Ward SOS expression Psos 01, the 
dipole free expression Bdf [4lj ]. or a non-perturbative 
finite difference approximation Bm p.[42| In the optimiza- 
tion code we use 0np- That is, we seek to maximize 
Pint = Bnp/Bmax- Once the optimization is complete, 
we use a comparison of 0np with /3sos and Bdf as a 
test of convergence. 

The exact computation of Bsos and Bdf requires sums 
over infinitely many states. We approximate them by 
summing over the 30 lowest energy levels. This is found 
to give an accurate value for the hyperpolarizability that 
is almost identical with the result we get when using up 
to 100 states. [H 

To illustrate the difference between the one- and two- 
electron cases, we take a closer look at the Schrodinger 
eigenvalue problems. For the one-electron problem in 
one-dimension, we solve 



h 2 d 2 * 



V(x)V = EV, 



(4) 



2m dx 2 

where the wave function ^ is a function of x. For two 
electrons without interactions, the Schrodinger equation 



s 2 * a 2 * 



2m \ dx\ 



dx\ 



+ {V(xx) + V(x 2 )) * = E9, (5) 



where $ is a function of x\ and x 2 . The solutions of 
Equation [5] are tensor products of solutions of Equa- 



tion |H and the energy levels are sums of energy levels 
of Equation @] The largest B values that we can obtain 
from Equation [5] are exactly twice those that we get from 
Equation[H but the largest value of Bint is lower by a fac- 
tor of exactly v2 because of the factor TV 3 / 2 that appears 
in Equation [2] for 0m AX- 

To include interactions between the electrons, we con- 
sider a model of the form 



— V 2 f + {V(x x ) + V(x 2 ) 
2m 

V 2 (\x 1 -x 2 \))^ = E^. 



(6) 



Here V(xi) is the external potential that acts on both 
electrons equivalently. It can be due to the nuclei or 
external influences such as an applied electromagnetic 
field. V^d^i — x 2 \) is the interaction potential between 
electrons, which is a sum of two components 

V 2 (\X! - x 2 \) = V 21 (\X! -x 2 \) + VziQxi - x 2 \), 

where V 2 \ and V 22 denote electrostatic repulsion between 
the electrons and spin interaction terms, respectively. 
These terms are piecewise linear, which is appropriate 
for a one-dimensional model. 

Defining d = \x± — x 2 \, we assume that the interaction 
potentials take the form 



V 21 (d) = h 2 .max C x 



L-d 



-.0 



(7) 



where L is the length of the molecule (we use L = 25 for 
all calculations); and, 



V 22 (d) = S ■ max j C 2 ^-j-^,0 



(8) 



where L = C3L with C3 -C 1 to ensure that the spin 
interaction acts only at close range. Note that this is a 
one-dimensional coulomb potential that satisfies the gen- 
eral condition \/ 2 V{xi—x 2 ) cx 5(xi— x 2 ), where 6(xi—x 2 ) 
is the Dirac delta function. In one-dimension, V = 
S is given by 



S = 



3ft 2 
4 



(triplet) 
(singlet). 



(9) 



Note that the parameters that we obtain during numer- 
ical optimization may be unphysical in the sense the in- 
teraction energies that result can be much larger than 
is typically observed between electrons; but, we stress 
that our goal is to understand the extremes of what is 
consistent with the Schrodinger equation. As discussed 
later, the resulting universal properties are unchanged if 
the parameters are constrained to a physically-reasonable 
range. 

The piece-wise potential energy function is the most 
general one-dimensional potential energy function inso- 
far as the pieces are small compared with the size of the 
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smallest features of the system. The form of the inter- 
action potential, on the other hand, is precisely the one- 
dimensional coulomb interaction. Thus, our model ac- 
counts for any system with interacting point charges, of 
which a two-electron one-dimensional molecule is a sub- 
set. 

For most of the experiments reported here we use 
Ci = 0.1, C-2 = 0.5, C3 = 0.1. These are ad hoc choices. 
We found that our results are not too sensitive to the 
choice of constants. However, if V21 is made too large 
relative to V22, the electrons are separated into two sep- 
arate wells and we get small values of Pint- In some 
experiments reported below, we allowed the constants to 
vary. We found that the spin interaction term V22 is 
much more important for making the intrinsic hyperpo- 
larizability larger than is the electrostatic repulsion term 
V 2l . 

Subspaces of symmetric and antisymmetric func- 
tions are invariant under any Hamiltonian of the form 
Given by Equation [6] Therefore we can consider 
the symmetric and antisymmetric eigenvalue problems 
separately. For the symmetric spatial wavefunctions 
(^>(x2,xi) = ^(xi, £2)) the spin part is a singlet, while 
for the antisymmetric spatial wavefunctions (^f(x2, x\) = 
— ty(x\, X2)) the spin part is a triplet. 

We approximate Equation [5] on a 20 x 20 grid of 
quadratic serendipity finite elements [43j . On the sym- 
metric subspace the number of degrees of freedom is cut 
nearly in half by the restriction ty(x2,xi) = ^(x\, x%). 
On the antisymmetric subspace it is more than halved. 
On each of the subspaces we use the implicitly restarted 
Arnoldi method [SiJ to compute the wave functions and 
energy levels. 

In all of our computations, the ground state is a sym- 
metric function (spin singlet), so the value of /3 depends 
only on the states that are symmetric in space. To see 
this, we consider the definition of the transition moments 
io m , which appear as factors in every term of the SOS 
or dipole-free expression for f3, 

XQrn — I I (tl+t 2 )^o(tl,t 2 )^ m (t 1 ,t2)dtldt2. 

Jo Jo 

If ^0 is symmetric and <P m is antisymmetric, then the 



integrand is antisymmetric and the integral xo m is zero. 
Thus this term contributes nothing to /3. 

In addition to the computations of Pint, we compute 
the matrix 32L \3&i 



'rap ~ u ra,p 



1 N 

n=0 



(Enm ~t~ E n p) X m n X n p, (1^) 



where N is the number of states used. Each matrix ele- 
ment of t( n \ indexed by m and p, is a measure of how 
well the (m,p) sum rule is obeyed when truncated to N 

states. If the sum rules are exactly obeyed, Tmp = for 
all m and p. We note that if the sum rules are truncated 
to an N-state model, the sum rules indexed by a large 
value of m or p (i.e. m,p ~ N) may be disobeyed even 
when the position matrix elements and energies are ex- 
act. We have found that the values of TmJ are small for 
exact wavefunctions when m < N/2 and p < N/2. So, 
when evaluating the r matrix to test our calculations, we 

consider only the components t^)n/2 p <n/2- 

Since the hyperpolarizability depends critically on the 
transition dipole moment from the ground state to the 
excited states, we use the value of Tqq ^ as one impor- 
tant test of the accuracy of the calculated wavefunctions. 
Additionally, we use the standard deviation of t^ n \ 



At (n) = 



v JV/2 ^N/2 / 
Z^m=0 Z^p=0 \^ 



rap 



N/2 



(11) 



which quantifies, on average, how well the sum rules are 
obeyed in aggregate, making a broader test of the 

accuracy of a large set of wavefunctions. 

We mentioned earlier that we did our computations on 
a 20 x 20 grid of finite elements. As a check on whether 
this mesh is fine enough, we repeated our final computa- 
tions on a 100 x 100 grid and found that we got the same 
results. 

Our code is written in MATLAB. 



III. RESULTS AND DISCUSSIONS 



We use various analytical functions as starting poten- 
tials, and then apply the numerical optimization tech- 
nique as described above to determine the optimal po- 
tential function for each. The tables that follow give the 
hyperpolarizability of the starting potential, (is followed 
by the hyperpolarizabilities of the optimized potentials. 
So, for example, when V = 0, we see that /3s = 0; but, 
the optimized potentials are all nonzero. 

Because the Nelder-Mead algorithm delivers only a lo- 
cal maximum, the final optimized potential can vary, de- 



pending on the choice of the starting potential. Table U 
summarizes the result of optimizing f3i n t for two inter- 
acting electrons with C\ — 0.1, C 2 — 0.5, and C3 = 0.1 
when starting from six different potentials, which are the 
same as in previous work 33] except that they have been 
re-scaled for convenience. Each starting potential (except 
V(x) = 0) was scaled so that V(L) = ?i 2 /10. Also, when 
we state that our starting potential is V(x) = x 2 , for ex- 
ample, we mean that it is the piecewise linear interpolant 
of x 2 . 
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TABLE II. Summary of results for the hyperpolarizability computed from the subspace of antisymmetric functions (triplet 
states). j3 a is the intrinsic hyperpolarizability of the starting potential while /3; nt is the value after optimization. The transition 
moments and energies are in the same dimensionless units as in Table [I] 

Function Js Prnt r^ 0) Ar (30) E 10 E 20 E 10 /E 20 x 00 x 10 JJfe 



V{x) 



0.6692 0.0014 0.0546 0.0528 0.1153 0.4582 5.7865 3.4477 0.7925 

x 0.0636 0.6690 0.0056 0.1472 0.0813 0.1817 0.4472 1.5841 2.7723 0.7902 

x 2 0.0057 0.6670 0.0025 0.1181 0.0650 0.1426 0.4556 2.4667 3.0848 0.7863 

y/x 0.1047 0.6682 0.0001 0.0999 0.0713 0.1582 0.4507 2.0296 2.9681 0.7926 

x + cos(x) 0.2113 0.6676 0.0023 0.0614 0.0782 0.1664 0.4698 4.4549 2.8294 0.7912 

tanh(x) 0.6691 0.0032 0.0570 0.0538 0.1174 0.4586 7.4199 3.4047 0.7900 



TABLE I. Summary of results for two interacting electrons - with Ci = 0.1, C 2 = 0.5, and C3 = 0.1 - using different starting 
potentials. /3 S is the intrinsic hyperpolarizability of the starting potential while fiint is the value after optimization. The 
transition moments and energies are in dimensionless units. To convert to specific units, consider an example where x nm is 
interpreted to be in units of angstroms. The energies would then be determined by multiplying all values of E n o by h 2 /ma 2 , 
with a = 10~ 10 m (1 A). In this case, the energy is in units of 1.2 x 10~ 18 J or about 7.6eV. 

Function Js p mt t (80) Ar (80) E w E 20 E 10 /E 20 x 00 x 10 

YM 

0.7068 0.0027 0.0704 0.0235 0.0488 0.4816 11.2267 -5.1493 0.7888 

x 0.6319 0.7063 0.0044 0.0574 0.0253 0.0528 0.4787 -11.8570 -4.9505 0.7872 

x 2 0.5523 0.7064 0.0035 0.0589 0.0249 0.0520 0.4799 -11.7233 -4.9913 0.7883 

y/x 0.6648 0.7063 0.0042 0.0599 0.0259 0.0540 0.4808 -12.0313 -4.8971 0.7888 

x + cos(x) 0.0311 0.7021 0.0044 0.0631 0.0381 0.0792 0.4808 -11.0082 -4.0336 0.7869 

tanh(a;) 0.0018 0.7068 0.0024 0.0601 0.0234 0.0488 0.4800 -11.2016 -5.1530 0.7886 

random 0.1066 0.7022 0.0080 0.0770 0.0395 0.0832 0.4749 14.6972 3.9642 0.7882 



It is important to stress that our approach is based on 
invariance under simple scaling. fj~3j so that the magni- 
tudes of the energies, dipole matrix elements, and hyper- 
polarizabilities are unimportant to the analysis. Rather, 
ratios of the observables, such as -E10/-E20, ^lo/^uf^ and 
P/Pmax are the key intensive quantities. As such, any row 
of Table Q] can be transformed by simple scaling to arbi- 
trarily adjust all the quantities, but leaving the intensive 
parameters unaffected. It is this property of simple scal- 
ing that makes possible a comparison of a broad range 
of molecular shapes and sizes using the intrinsic hyper- 
polarizability as the metric. 

The values of Pi n t reported in the table were computed 
using font = Pnp I 'ft max- When we checked the com- 
putation using (3 sos in place of (3 NP , we got the same 
values to at least four decimal places, thus, we choose to 
use Pnp in all of our optimization calculations. When 
we used Pd f , we had agreement to within a fraction of a 
percent. In the last column of Table [H x^ ax is the mag- 



nitude of the fundamental limit of the position matrix 
element X10 for a two-electron system, and is given by, 

max _ ti /if)N 
10 _ / — £5 — V I 

\Jmkiva 

As was found for the one-electron case,[l3|, H3| a quan- 
tum system with a hyperpolarizability near the funda- 
mental limit shares certain universal properties. For ex- 
ample, all optimized values of Pint are just over 0.7, for 
which the energy ratio -Bio/^o is just under 0.5 and the 
intrinsic transition moment is about 0.79. Note that 
if electron interaction as characterized by the C parame- 
ters is decreased smoothly to zero, the optimized hyper- 
polarizability approaches 0.5, as is predicted analytically 
for the non-interacting limit. When electron interactions 
are made very large, the optimized intrinsic hyperpolar- 
izability does not increase beyond Pint = 0.708. 



While we have focused mainly on well-defined potential 
energy functions, it is also possible to apply Monte Carlo 
methods to randomly pick the shape of the potential. In 



the first simulation, we picked the random numbers to be 
normally distributed with mean zero and standard devi- 
ation 1. The optimization algorithm got trapped in the 
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parameter space where Pint was small. However, for a set 
of random numbers with the standard deviation reduced 
to 0.1, the simulation gave Pint ~ 0.7. The resulting 
parameters are shown in the last row of Table fl] Once 
again, the same universal parameters are observed as for 
the other optimized potentials. 

It is not likely that the universal properties that 
we observe are coincidental or artifacts of the calcu- 
lation given independent calculations using widely dif- 
ferent approaches. In addition to the large number of 
one-dimensional potentials for one and two electrons, 
for which these same universal properties result. [32l. 
[33[ the same universal properties are found when al- 
lowing point nuclei to move in a plane to form a 
molecule, [45[ applying electromagnetic fields to two- 
dimensional molecules, [37j and using constrained po- 
tentials with properties more similar to real molecules 
and including continuum states. 146] Monte Carlo calcu- 
lations, which probe a much larger parameter space than 
the former examples are also observed to share some uni- 
versal properties, though they are less restrictive than in 
the cases where a potential energy function and vector 
potential is used. 

In all simulations, the energies of the triplet states were 
tracked and found to be of higher energy than the ground 
state singlet state. Since the triplet states are orthogonal 
to the space of singlet states, no triplet states contribute 
to the results. To investigate the hyperpolarizability due 
to the triplet states, we used the lowest energy triplet 
state to calculate the dipole moment from which the hy- 
perpolarizability is determined using the nonperturba- 
tive method. We also used the SOS expressions with the 
lowest energy triplet state as the ground state to cal- 
culate the hyperpolarizability. Thus, these calculations 
only consider transitions between triplet states. Table [TT1 
summarizes the results. 

In all cases, the optimized intrinsic hyperpolarizability 
for the triplet subspace is just short of 0.67, which is lower 
than the value of 0.71 found in the singlet subspace. Sim- 
ilarly, the other parameters, such as the energy ratio and 
normalized transition moment, are close but not identical 
to the singlet values. Thus, we conclude that the qualities 
of the optimized system using only triplet states is simi- 
lar to those of the singlet subspace. However, when using 
the true ground state of the system, the same universal 
properties are observed for two interacting electrons as is 
observed for a single electron. 

We define j3 nm as the contribution from states n and 
m to the SOS sum,[33| 



trix elements are at least 10 times smaller. Thus, the non- 



Pn 



(13) 



„ iTn . The hyperpolarizability is 
,, m over all excited states of the 



where x nm = x nm ~ XooS. 
calculated by summing /3 
system. Figure [T] shows a plot of P nm that is representa 
tive for all the potentials plotted in Figure |U Typical val- 
ues are Pn = —0.9, P21 = P12 — 0.09, and all other ma- 




FIG. 1. Surface plot of j3 nrn that is representative of all the 
optimized hyperpolarizabilities in Table [I] 



interacting electron potentials that optimize pint lead to 
hyperpolarizabilities that obey the three-level ansatz, [32j 
that is, a system with /3 near the fundamental limit is de- 
scribed by a three-level model. 

The three-level ansatz has been shown to be univer- 
sally obeyed in systems described by a one-dimensional 
potential, [33 point nuclei in a plan e. l4o| and for one- 
dimensional monte-carlo simulations. 38] In the case of 
arbitrarily applied electric fields to point nuclei in two- 
dimensions, the optimized diagonal component of the hy- 
perpolarizability obeys the three-level ansatz. However, 
optimized p xxy yields p int = Pxxy/P™™ » 0.9 and the 
intrinsic transition moment is 0.51. In this case, the hy- 
perpolarizability is dominated by two excited states; but, 
an additional pair of excited states contribute about 25%. 

In further experiments, starting with the potential de- 
picted in Figure [2J we attempted to increase Pint even 
further by varying the constants Ci, C2, and C3 appear- 
ing in the interaction terms (Equation [7] and Equation 
[5]). These constants are all constrained to be positive. 
We achieved only a modest increase from (3 int = 0.7068 
to j3i n t — 0.7080, still no better than the largest intrin- 
sic hyperpolarizability attained by varying the potential 
energy function; [33] but the largest found in our simula- 
tions. The final values of Ci, C2, and C3 were 0.0034, 
0.8787, and 0.4942. The strength of the electrostatic re- 
pulsion term is governed by C\, which decreased from 
0.1 to 0.0034 during the optimization. Similar results 
were seen in other experiments of this type. Thus the 
electrostatic repulsion term seems not to be important 
for maximizing Pint- The strength of the spin interac- 
tion term is governed by C2, which increased from 0.5 to 
0.8787 during optimization. This suggests that the spin 
interaction term is important for maximizing Pint- 
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TABLE III. Parameters resulting from optimizing the intrinsic hyperpolarizability starting with the potentials shown and 
varying the strengths of the electron interaction terms. For the V(x) = starting potential, both the electron interaction 
parameters and potential function were varied. f5 s is the intrinsic hyperpolarizability of the starting potential while fimt is the 
value after optimization. 



Function /3 S Pint r ( * 0) Ar (80) E 10 E 20 
V(x) 


E10/E20 


XlO 


2:10 


Figured 0.7068 0.7080 0.0018 0.0627 0.0235 0.0490 
x 0.6319 0.6663 0.0069 0.0149 0.0442 0.0804 
0.7074 0.0010 0.0549 0.0249 0.0516 


0.4803 
0.5494 
0.4827 


5.1265 
-4.1252 
-4.9887 


0.7865 
0.8670 
0.7870 



TABLE IV. Electron interaction parameters resulting from 
optimizing the intrinsic hyperpolarizability starting with the 
potentials shown and varying the strengths of the electron in- 
teraction terms. For the V(x) = starting potential, both the 
electron interaction parameters and potential function were 
varied. p s is the intrinsic hyperpolarizability of the starting 
potential while Pint is the value after optimization. 



Function p s Pint 


Ci 


C-2 C3 


V(x) 






Figured 0.7068 0.7080 


0.0034 


0.8787 0.4942 


x 0.6319 0.6663 


0.000087 


14.75 0.0612 


0.7074 


0.0933 


0.5283 0.1033 



In another experiment we used the potential V(x) = x 
and started with C\ = 0.1, Ci = 0.5, and C3 = 0.1. 
Holding V(x) fixed, we tried to increase f3i nt by varying 
Ci, C2, and C3. With this potential we already have a 
rather good @i n t = 0.6319 to begin with. By varying the 
constants we were only able to realize a modest in- 
crease in I3i n t to about 0.6663, which was obtained with 
Ci = 0.000087, C 2 = 14.75, and C 3 = 0.0612. With such 
a small C\ , the electrostatic repulsion term is essentially 
zero. The large value of Ci means that the spin interac- 
tion term is very strong, and the small value of C3 means 
that the spin interaction operates only over a very short 
range. We conclude that varying the electron interaction 
potential alone does not have a large effect on improving 
the intrinsic hyperpolarizability. 

In another experiment we started with the potential 
V{x) = and the initial constants C\ = 0.1, Ci = 0.5, 
and C3 =0.1. When allowing both the potential energy 
function and electron interaction parameters to vary, we 
obtained f3i n t — 0.7074. The final values of the electron 
interaction parameters were C\ — 0.0933, Ci = 0.5283, 
and C3 = 0.1033, so they changed little from their initial 
values. 

Table Hill summarizes some of the parameters found for 
the optimized intrinsic hyperpolarizabilities when elec- 
tron interactions are present. In the two cases where 
the intrinsic hyperpolarizability is greater than 0.7, the 
same universal values are observed; i.e. the energy ratios 
and the intrinsic transition moments are almost identical 
to the ones in Table HI Furthermore, (3 nm for all three 



optimized potentials are found to be dominated by only 
one pair of excited states. Thus, the three-level ansatz is 
again obeyed when electron interactions are important. 

Based on these observations, we conclude that when 
the diagonal component of the hyperpolarizability is op- 
timized, all systems - including ones with two interacting 
electrons - share the same universal properties; and, the 
three-level ansatz holds. This observation suggests that 
when investigating the properties of a quantum system, 
the numerical optimization of the simplest models can be 
used to identify universal properties. This paradigm is 
in stark contrast to studies that seek to understand the 
properties of a specific system, pi |47]| and is more akin to 
the theoretical studies by Meyers et al, who studied the 
effect of an electric field on polyenes [48| which confirmed 
the concept of bond length alternation. [4^, [5(| However, 
our approach of optimizing the nonlinear response using 
the potential energy function is even more general. The 
only constraint we apply is that the sum rules are obeyed, 
which demands consistency with quantum mechanics but 
is otherwise general. 

Arguably, the most important result of universal prop- 
erty studies is that there appear to be a very large num- 
ber of ways to make a system with a nonlinear-optical 
response near the fundamental limit with Pi nt = O.T. [l~3| 
Szafruga and coworkers performed numerical optimiza- 
tion studies under the constraint that the potential en- 
ergy varies by an amount that is typical in a molecule, 
as well as extending the computational space beyond the 
potential well to approximate continuum states. [46j It 
was found that with those restrictions the same universal 
properties were observed, and (3i nt ~ 0.708. Therefore, it 
appears that it should be possible to make many classes 
of yet-to-be-found molecules that are at the fundamental 
limit. 

Some of the parameters shown in Table IIVI yield elec- 
tron interaction energies that are larger than typically 
observed. However, even in the strong interaction limit, 
the effect on the optimized hyperpolarizability is not 
large. Thus, we conclude that when electron interac- 
tions are included, the universal properties remain the 
same and that electron interactions may often be unim- 
portant when studying universal properties of optimized 
systems. This is not to say that electron interactions are 
unimportant when calculating the properties of a specific 
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FIG. 2. Optimized potential energy function and electron 
density for the lowest-energy 15 states of a 30-state model. 
Starting potential is V(x) = 0. 



molecule. When we limited the parameter space in the 
present studies in the way described in Szafruga[46j|. the 
observed universal properties with electron interactions 
remained the same. 

Finally, we compare the charge densities that result for 
the optimized potential energy function with two elec- 
trons with prior results with one electron. Figure [2] and 
Figure [3] show the optimized potential energy functions 
for starting potentials V(x) = and V(x) = tanh(a;). 
Also shown is the electron density for the lowest-energy 
fifteen states of a thirty -state model. Note that the nu- 
merical labels on the y axes refer only to the potential 
energy functions and all electron densities are positive. 
This is a common observation that very different poten- 
tial energy functions that lead to very different wavefunc- 
tions always share the same universal properties. Thus, 
there is a unifying underlying structure that is shared 
by these optimized systems that are not apparent simply 
from the electron densities. 

In past studies. [33ll42j it was found that two classes of 
wavefunction lead to optimized hyperpolarizability: (1) 
Localized charge density resulting from randomly oscil- 
lating potentials, reminiscent of Anderson Localization 
for the tanh(x) starting potential, and (2) delocalized 
wavefunctions. In both cases, two excited states domi- 
nated in their contribution to the hyperpolarizability. In 
the present studies, there is a hint of localization in the 
ground and first excited states, but the others are delocal- 
ized. However, this result may be due to the fact that the 
present calculation used piecewise continuous potentials 
rather than splines, which were used in previous studies. 

To study localization, it is more instructive to consider 
the full two-electron wavefunctions. Figure [4] shows the 
electron density \ipn{ x i, x 2)\ 2 as a function of X\ and x^ 
for an optimized quantum system using the starting po- 
tential V(x) = tanh(x). The electrons are correlated in 



FIG. 3. Optimized potential energy function and electron 
density for 15 states. Starting potential is V(x) — tanh(x). 
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FIG. 4. Contour plots of electron density functions for the 
nine lowest energy levels (left to right, top to bottom) for the 
optimized potential with starting potential V = tanha;. 



the four lowest energy eigenstates. Many of the higher 
states show that the two electrons are localized in dif- 
ferent spatial regions. This is reminiscent of the one- 
dimensional one-electron potential, where the wave func- 
tions do not overlap. Note that the electron densities 
appear similar in the excited states of all optimized sys- 
tems. 

The transition dipole moment amplitude from the 
ground state, tp n {xi,X2) {xi + x%) ipo(xi,X2), is shown in 
Figure [S] The only transition dipole moment densities 
that are large are for transitions from the ground state to 
the first and second excited states. For the higher states, 
the transition moment amplitudes are either small; or, 
large positive regions are balanced by small negative 
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FIG. 5. Contour plots of transition dipole moment ampli- 
tude tp n (xi, X2) (xi + x%) ipo{x\, X2) with n = 1 to 9 (left to 
right and top to bottom) for the optimized potential energy 
function with starting potential V = tanha;. 



regions. Thus, the three-level ansatz is obeyed either 
through minimal overlap between the different eigenfunc- 
tions, or cancelation between positive and negative prob- 
ability amplitude. 

One might ask how the potential energy function can 
be manipulated to make a real system with an opti- 
mized nonlinear optical response. There are many ar- 
tificial structures, such as nanoparticles, quantum wires, 
and multiple quantum wells in which there is some con- 
trol over the structure and the resulting potential profile. 
For example, nanowires can be used as building blocks for 
nanoscale electronic and optoelectronic devices, [5l| and 
multiple quantum wells can be made, for example, with a 
band edge that can be used to make a s yste m with a large 
quadratic electroabsorption coefficient. [52j Recently, Hall 
and coworkers used graded barriers to control the optical 
properties of ZnO/ZnMgO quantum wells with an in- 
trinsic internal electric field. [531] Clearly, there are many 
alternatives to molecules that can be used to deliberately 
engineer a potential profile. 

Wang and coworkers have shown that the potential 
energy profile of a molecule can be approximated by us- 
ing a linear combination of atomic potentials. [54] Thus, 
the results of our work may be able to be connected 
to real molecules using such an approach. We stress 
that optimization of the potential energy function is 
a more general approach than other numerical opti- 
mization techniques, with the exception of the Monte 
Carlo techniques, [H, HI] because potentials apply to the 
broader set of systems of which molecules are a small 
subset. Our present work shows that the observed uni- 
versal properties are not restricted to just one electron 
systems. 

It is no simple matter to determine the nature of the in- 



teraction between electrons that leads to the iV 3 / 2 scaling 
law. The fact that the fundamental limit of the hyper- 
polarizability scales as N 3 / 2 follows from the sum rules 
and the limits they set for the oscillator strengths. If the 
full oscillator strength resides in one state, the modulus 
squared of the transition moment to that state from the 
ground state is proportional to the number of electrons, 
thus, the transition moments scale as the square root 
of the number of electrons. Since the hyperpolarizabil- 
ity scales as the third power of the transition moment, 
it scales as TV 3 / 2 . Note that the sum over the dipole 
terms of the form Ax n \x n o\ 2 can be re-expressed using 
the sum rules in terms of sums of the products of the 
form xon^nm^mO-IHl Thus, while counterintuitive, the 
dipole terms acting together will also scale as iV 3 / 2 . 

In the absence of interactions, the optical properties 
are obviously additive over the molecular units. When 
allowed to interact, the largest power law allowed by the 
fundamental limits is TV 3 / 2 , and the optimization algo- 
rithm finds this local maximum. The fact that the pa- 
rameters that characterize the interaction show no uni- 
versal behavior (see Table IIV[) implies that there is no 
unique interaction that leads to the optimal response. 
Thus, the question about the nature of the interactions 
that leads to optimized scaling does not have a unique 
answer. 

In closing, we note that the focus of our present work is 
on the off-resonant hyperpolarizability. Since the funda- 
mental limit of the dispersion of the hyperpolarizability 
has also been calculated, (5|| it would be a straightfor- 
ward matter to calculate the potential energy function 
that is optimized at any set of wavelengths. Such calcula- 
tions would be useful in determining the optical material 
properties required for a particular device. Since there 
are an infinite number of combinations of wavelengths, 
we have limited our present studies to the off-resonance 
regime. 

IV. CONCLUSIONS 

Rather than calculating the nonlinear-optical response 
of large numbers of complex molecules and searching for 
correlations between structure and nonlinear response, 
our work focuses on developing general principles that 
broadly apply to all quantum systems; and, endeavors 
to identify universal properties that can be applied to 
build an understanding of what leads to a large nonlinear- 
optical susceptibility. 

The present work seeks to determine whether or not 
electron interactions and Fermi correlations change the 
universal behavior of an optimized system as is observed 
for one electron systems. We have found that the univer- 
sal behavior is indeed the same. The intrinsic transition 
moment, energy ratio of the first two dominant excited 
states, and the intrinsic hyperpolarizabilities for an opti- 
mized system are all the same for a broad range of start- 
ing potentials. Also, the three-level ansatz is again found 
to be obeyed. Since electron interactions appear to have 
a small effect on the conclusions derived from numerical 
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optimization, simple one-electron models may be enough 
to study the broad principles underlying light-matter in- 
teractions. 

Optimization studies are inherently tricky because one 
can never be certain whether or not a global maximum 
has been found. However, as we continue to gener- 
ate more numerical data from a broader range of start- 
ing potentials, molecular symmetries, externally-applied 
electromagnetic fields, and electron interactions, we con- 
tinue to observe the same universal properties - leading 
to the tentative conclusion that perhaps the data as a 



whole is representative of the properties of the true global 
maxima. Similarly, with only a handful of exceptions, 
experimental work that focuses on the development of 
new molecules finds that the hyperpolarizabilities are far 
short of the fundamental limit. The broad principles that 
we seek to develop aims to bridge the gap between the 
limits and real molecular systems. Clearly, altogether 
new paradigms will need to be developed to attain what 
is possible. 
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